#define INF     2e10f
#define DIMX     600
#define DIMY     300
#define dx      (4.0f/DIMX)
#define dy      (2.0f/DIMY)
#define X0  -2.0f
#define Y0  -1.0f
#define Z0  -1.0f
#define ox 0.0f
#define oy 0.0f
#define oz 0.0f 
#define ENTITYS 2

#define ABS(x) (x>0.0f ? x:-x)
#ifndef MAX
#define MAX(a,b) (a > b ? a : b)
#endif
#ifndef MIN
#define MIN(a,b) (a > b ? b : a)
#endif

// //////////////////产生随机数///////////////////////////////////////////////////
/*
Part of MWC64X by David Thomas, dt10@imperial.ac.uk
This is provided under BSD, full license is with the main package.
See http://www.doc.ic.ac.uk/~dt10/research
*/
#ifndef dt10_mwc64x_rng_cl
#define dt10_mwc64x_rng_cl

//! Represents the state of a particular generator
typedef struct{ uint x; uint c; } mwc64x_state_t;

enum{ MWC64X_A = 4294883355U };
enum{ MWC64X_M = 18446383549859758079UL };

void MWC64X_Step(mwc64x_state_t *s)
{
	uint X=s->x, C=s->c;
	
	uint Xn=MWC64X_A*X+C;
	uint carry=(uint)(Xn<C);				// The (Xn<C) will be zero or one for scalar
	uint Cn=mad_hi(MWC64X_A,X,carry);  
	
	s->x=Xn;
	s->c=Cn;
}
//! Return a 32-bit integer in the range [0..2^32)
uint MWC64X_NextUint(mwc64x_state_t *s)
{
	uint res=s->x ^ s->c;
	MWC64X_Step(s);
	return res;
}

#endif
inline float intersect_plane(float3 O, float3 D, float3 P,float3 N)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // plane (P, N), or +inf if there is no intersection.
    // O and P are 3D points, D and N (normal) are normalized vectors.
    float denom = dot(D, N);
    if (ABS(denom) < 1e-6)
        return INF;
    float d = dot(P - O, N) / denom;
    if (d < 0)
        return INF;
    return d;
}

inline float intersect_sphere(float3 O,float3 D,float3 S,float R)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // sphere (S, R), or +inf if there is no intersection.
    // O and S are 3D points, D (direction) is a normalized vector, R is a scalar.
	float3 SO=S-O;
    float b = dot(D, SO);
    float c = dot(SO, SO) - R * R;
    float disc = b * b - c;
    if (disc > 0 && b>0){
		float rt=sqrt(disc);
        return ((b-rt)>0 && (b+rt)>0)?(b-rt):(b+rt); //同时处理反射和折射
	}
    return INF;
}

inline float intersect_sphere2(float3 O,float3 D,float3 S,float R)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // sphere (S, R), or +inf if there is no intersection.
    // O and S are 3D points, D (direction) is a normalized vector, R is a scalar.
	float a = dot(D, D);
    float3 OS = O - S;
    float b = 2 * dot(D, OS);
    float c = dot(OS, OS) - R * R;
    float disc = b * b - 4 * a * c;
    if (disc > 0){
        float distSqrt = sqrt(disc);
        float q =   b < 0 ? (-b - distSqrt) / 2.0 : (-b + distSqrt) / 2.0;
        float t0 = q / a;
        float t1 = c / q;
        t0= MIN(t0, t1);
		t1 = MAX(t0, t1);
        if (t1 >= 0)
            return t0 < 0 ? t1 : t0;
	}
    return INF;
}

inline float intersect(float3 O,float3 D,float4 obj,float3 normal,__global float4* debugarr)
{
	//半径无穷大，平面
    if (obj.w == INF)
	{
		// debugarr[0]=(float4)(INF,INF,INF,INF);
		//debugarr[0].xyz=normal;
		float t=intersect_plane(O, D, obj.xyz, normal);
		//debugarr[0].w=t;
		//debugarr[1].xyz=O;
		//debugarr[2].xyz=D;
        return intersect_plane(O, D, obj.xyz, normal);
	}
    else
        return intersect_sphere(O, D, obj.xyz, obj.w);
}

inline float find_intersect_object(float3 O,float3 D , int* obj_idx,__global const float4* pos,__global const float4* normal,__global float4* debugarr)
{
    float t = INF;
    *obj_idx = -1;
    for (int i=0;i< ENTITYS;i++){
		//debugarr[i]=normal[i];
		// debugarr[i].w=i+1;
        float t_obj = intersect(O, D, pos[i],normal[i].xyz,debugarr);
        if (t_obj < t){
            t= t_obj;
			*obj_idx =i;
		}
	}
    return t;
}

inline int find_shadow(float3 M,float3 N,float3 toL,int obj_idx,__global const float4* pos,__global const float4* normal,__global float4* debugarr)
{
	float ll=length(toL);
    for(int k=0;k< ENTITYS;k++){
        if (k != obj_idx){
            float t = intersect(M + N * 0.0001f, normalize(toL), pos[k],normal[k].xyz,debugarr);
            if(t<=ll) return 1;
		}
	}
    return 0;
}

inline float3 get_normal(int obj,float3 M ,__global const float4* pos,__global const float4* normal)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return normal[obj].xyz;
    else{
		return normalize(M-p.xyz);
	}
}

inline float3 get_color(int obj,float3 M ,__global const float4* pos,__global const float4* color)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return (float3)(((int)(M.x*2+20)%2)==((int)(M.z*2+20)%2) ? 1.0f : 0.0f);
    else{
		return color[obj].xyz;
	}
}

//计算反射光线
inline float3 reflect(float3 v, float3 n)
{
	return v - 2 * dot(v, n) * n;
}

// 计算折射光线
bool refract(float3 v, float3 n, float ni_over_nt, float3 refracted)
{
	float3 uv = normalize(v);
	float dt = dot(uv, n);
	float discriminant = 1.0 - ni_over_nt * ni_over_nt * (1 - dt * dt);
	if (discriminant > 0)
	{
		refracted = ni_over_nt * (uv - n * dt) - n * sqrt(discriminant);
		return true;
	}
	else
		return false;
}

float3 random_in_unit_sphere(mwc64x_state_t *rng) {
    float3 p;
    do {
        p = 2.0f*(float3)((MWC64X_NextUint(rng)/4294967296.0f), (MWC64X_NextUint(rng)/4294967296.0f), (MWC64X_NextUint(rng)/4294967296.0f)) - (float3)(1.0f,1.0f,1.0f);
    } while (length(p) >= 1.0f);
    return p;
}

__kernel void ray_tracing( __global const float4* color,__global const float4* pos,
						   __global const float4* normal,__global const float4* material,
						   __global const float4* lights, __global unsigned char * bitmap,
						   __global float4* debugarr
						   )
{
    // map from threadIdx/BlockIdx to pixel position
    int i = get_global_id(0);
    int j = get_global_id(1);
	if (i>=DIMX || j>=DIMY) return;
    int offset = i + j * get_global_size(0);
	// 随机采样 抗锯齿
	mwc64x_state_t rng={i*9000+19797,j*1040+1432};
	MWC64X_NextUint(&rng);
    float3 O =(float3)(ox,oy,oz); //眼睛的位置
	
	
	int obj_idx=0;// 与视线相交的实体的编号
	// 一条视线最多追踪4次
	float3 color_ij=(float3)(0,0,0);
	
	int ns=100;//采样次数
	//float3 tmp=(float3)(0.0f);
	for (int s=0;s<ns;s++){
		float reflection=1.0f;
		float3 tmp=(float3)(0.0f);
		float3 Q =(float3)(X0+((MWC64X_NextUint(&rng)/4294967296.0f)+i)*dx,Y0+(MWC64X_NextUint(&rng)/4294967296.0f+j)*dy,Z0); //视线和屏幕的交点
		//float3 Q =(float3)(X0+(i)*dx,Y0+(j)*dy,Z0); //视线和屏幕的交点
		float3 D=normalize(Q-O);// 视线的方向
		float3 ORG=O;
		for (int n=0;n<100;n++){
			float dt= find_intersect_object(ORG, D,&obj_idx,pos,normal,debugarr);
			if( dt == INF) {// 不相交直接退出
			float t = 0.5f*(D.y + 1.0f);
			tmp=((1.0f-t)*(float3)(1.0f, 1.0f, 1.0f) + t*(float3)(0.5f, 0.7f, 1.0f))*reflection;
			break;
			}else{
			// 交点
			float3 M = ORG + D * dt;
			//debugarr[3].xyz=M;
			// 法线
			float3 N = get_normal(obj_idx,M ,pos,normal);
			// // 指向光源
			// float3 L=lights[0].xyz;
			// float3 toL =L - M;
			// // 指向眼睛
			// float3 toO = normalize(O - M);
			// // 判断光源是否被挡住
			// int shadowed=find_shadow(M,N,toL,obj_idx,pos,normal,debugarr);
			// toL=normalize(toL);
			// // 材质参数
			// float4 mat=material[obj_idx];
			// float3 temp=(float3)(0,0,0);
			// // 计算环境光
			// temp+=0.05f;
			// if (!shadowed){
			// // 计算漫反射 Lambert shading (diffuse).
				// temp+=mat.y*MAX(dot(N, toL), 0) * get_color(obj_idx,M ,pos,color);
			// // 计算镜面反射 Blinn-Phong shading (specular).
				// temp+=mat.z*pow(MAX(dot(N, normalize(toL + toO)), 0), 50) * lights[1].xyz;
			// }
			// 按反射率计算贡献值
			//tmp=(N+1.0f)/2.0f;
			reflection*=0.5f;
			// 计算将要追踪的下一条反射光线
			D=normalize(N+random_in_unit_sphere(&rng));
			ORG=M+0.0001f*N; //反射
			}
		}
		color_ij+=tmp;
	}
	color_ij/=ns;
    bitmap[offset*4 + 0] = (int)(MIN(color_ij.x,1) * 255);
    bitmap[offset*4 + 1] = (int)(MIN(color_ij.y,1) * 255);
    bitmap[offset*4 + 2] = (int)(MIN(color_ij.z,1)  * 255);
    bitmap[offset*4 + 3] = 255;
}